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ABSTRACT 

This Letter presents the first results from the observations of LS I +61° 303 using Large Area Telescope data 
from the Fermi Gamma-Ray Space Telescope between 2008 August and 2009 March. Our results indicate 
variability that is consistent with the binary period, with the emission being modulated at 26.6 ±0.5 days. This 
constitutes the first detection of orbital periodicity in high-energy gamma rays (20 MeV-100 GeV, HE). The 
light curve is characterized by a broad peak after periastron, as well as a smaller peak just before apastron. 
The spectrum is best represented by a power law with an exponential cutoff, yielding an overall flux above 100 
MeV of 0.82 ± 0.03(stat) ± 0.07(syst) 10" 6 ph cm" 2 s _1 , with a cutoff at 6.3 ± l.l(stat) ± 0.4(syst) GeV and 
photon index T = 2.21 ± 0.04(stat) ± 0.06(syst). There is no significant spectral change with orbital phase. 
The phase of maximum emission, close to periastron, hints at inverse Compton scattering as the main radiation 
mechanism. However, previous very high-energy gamma ray (>100 GeV, VHE) observations by MAGIC and 
VERITAS show peak emission close to apastron. This and the energy cutoff seen with Fermi suggest the link 
between HE and VHE gamma rays is nontrivial. 

Subject headings: binaries: close — stars: variables: other — gamma rays: observations — X-rays: binaries 
— X-rays: individual (LS I +61 °303) 
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1. INTRODUCTION 

The high-mass X-ray binary LS I +61°303 (=V615 Cas) 
has long been plausibly associated with a high-energy (HE, 
20 MeV-100 GeV) gamma-ray source, although never before 
confirmed. The disco very of the COS B source 2CG 135+01 
dHermsen et alj[l977h quickly brought attention to this binary 
system's Be star localized within its error box, becau se of its 
unusual periodic radio emission dGregorv et aill 19791) and its 



rego 

X-ray emission dBignami et al.l 1 1 981D . 2CG 135+01 was to 
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remain one of the brightest sources known in the HE gamma- 
ray sky, with a flux of ~ 10~ 6 ph s" 1 cm" 2 above 100 MeV 
dSwanenburg et al.lll981l) . In the 1990s, EGRET detected the 
source with high confidence at the same average flux level 
and derived a powe r-law photon index of T = 2.05 ± 0.06 
dKniffen et alHl997l) . Although there are no other objects of 
note (ra dio-loud AGN , or pu lsars) coinciding with the 3EG 
source (Hartma n et al.l [19991) . its positional uncertainty was 
not small enough to firmly associate the gamma-ray source 
with the binary. Variability in the EGRET light curve could 
be neither fir mly established nor related to variabil ity at other 
wavelengths dTavani et al.ll998tlNolan et al.l2003l) . Recently, 
AGILE has reporte d detecting the source at the same flux level 
dPittori et al.ll2009l) . 

LS I +61°303 is an unusual binary system exhibiting strong 
variable emission from the radio to X-ray and TeV energies. 
At radio wavelengths the source has been shown to exhibit 
radio outbursts that are modulated on an orbital period of 
26.49 60 ± 0.0028 days dTavlor & Gregory! [l98l iGregorvl 
|2002|). The phase of radio maximum has also been shown 
bv lGregorvl d2002l) to vary with a super-orbital period of 1667 
± 8 days. Observations of orbital modulation in the optical 
place constraints on the binary system parameters. The bi- 
nary has an eccentric orbit (e=0.55-0.72) and the Be star ra- 
dial velocity is consistent with a neutron star companion or, 
if the orbital inclination is < 25°, with a > 3M ffi black hole 
dHutchings & Cramptonlfl98Tl: ICasares et alj|2005l) . Signifi- 
cant uncertainty still exists in key parameters of the orbital 
soluti on of the system (Grundstrom et al. 2007; Aragon a"et al.l 
120091) . 

Behavior in the X-ray band is much more complicated. 
Orbital modulation has been reporte d with the peak of 
emission appearing at phases 0.6-0.7 (Par edes et al.l Il99~7l : 
Esposito et al. 2007). However the modulation is not 
smooth, with short timescale flares a nd very strong orbit- 
to-orbit variability (Smi th et all 12009). Broad-band spec- 
tral analysis of XMM-Newton and I NTEGRAL data by 
IChernvakova. Neronov, & Walter! d2006l) reveal LS I +6 1 303 
to be well fitted by a simple absorbed power-law with a hard 
photon index, T ~1.5, in the 0.5-100 keV band. 

The MAGIC telescope detected a variable very high- 
energy (VHE >100 GeV) gamma-ray source coincident 
with LS I +61°303 dAlbert et all 120061) : a result that has 
been independently con firmed by the VERITAS collabora- 
tion dAcciari et alj|2008l) . More recently, the MAGIC collab- 
oration has further reported that the VHE emis sion is peri- 
odic a t the 26.5 day orbital period of the system dAlbert et alj 
2009). The VHE emission is consistently highest close to 
apastron, when the compact object is farthest from the Be 
star, and remains und etected at perias t ron. L ike LS 5039 
and PSR B 1259-6 3 dAharonian etail l2005allbl) . and con- 
trary to Cyg X-l dAlbert et alj I2007D . LS I +61° 303 is a 
gamma-ray binary with its spectral energy distribution peak- 
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ing in HE gamma rays (for a full SE P see iDubusI l2006t 
IChernyakova. Neronov. & Walterjl2006l) . 

Calculations of the theoretical expectations of the gamma 
ray emissi on from LS I +61°303 go back almost three 
decades ( Mara schi & Trevesl 1 198 lh . and there has been 
a recent burst of activity following the MAGIC detec- 
tion. Two scenarios have been put forward involving ei- 
ther the relativistic wind of a young, rotation-powered pul - 
sar dDubusl2006HSierpowska-Bartosik & Torresll2008l 120091) . 
or the re lativistic jet of an accreting black ho l e or neu- 
tron star (iRomero. C hristianse n. & Orel lana 20051 Bednarekl 
l2006HGupta & Bottcherll2006t iBosch-Ramon et al.H2006l> . In 
light of the orbital modulations seen in radio, X-ray, and VHE 
gamma rays, a detailed light curve in the HE gamma ray do- 
main (where most of the energy is output) is an essential piece 
to identify the main radiative process at work and model the 
source. 

2. DATA REDUCTION AND RESULTS 

The Fermi Gamma-ray Space Telescope was launched on 
2008 June 11, from Cape Canaveral, Florida. The Large Area 
Telescope (LAT) is an electron-positron pair production tele- 
scope, featuring solid state silicon trackers and cesium iodide 
calori meters, sensitive to photons from ~ 20 MeV to > 300 
GeV (lAtwood et alJ 120091) . Relative to earlier gamma-ray 
missions the LAT has a large ~ 2.4 sr field of view, a large ef- 
fective area (~ 8000 cm 2 for > 1 GeV on axis) and improved 
angular resolution or point spread function (PSF, better than 
1° for 68% containment at 1 GeV). The Fermi survey mode 
operations began on 2008 August 4, after the conclusion of 
a flawless commissioning period. In this mode, the observa- 
tory is rocked north and south on alternate orbits to provide 
more uniform coverage so that every part of the sky is ob- 
served for ~30 minutes every 3 hours. Thus Fermi is ideally 
suited for long term all-sky observations. The dataset for this 
analysis spanned 2008 Aug 4, through 2009 Mar 24. Thus 
LS I +61° 303 was observed for approximately 9 orbital peri- 
ods. 

The data were reduced and analysed using the Fermi Sci- 
ence Tools v9r8 package 59 . The standard onboard filtering, 
event reconstruction, and classification were applied to the 
data (lAtwood et al.l l2009h . and for this analysis the high- 
quality ("diffuse") event class is used. Time periods when the 
region around LS I +61°303 was observed at a zenith an- 
gle greater than 105° were also excluded to avoid contami- 
nation from Earth albedo photons. With these cuts, a pho- 
ton count map of a 10° region around the binary is shown 
in Fig. Q] The alignment of the LAT pointing direction with 
the celestial frame was calibrated using a large set of high 
latitud e gamma-ray sources to better than 10" dAbdo et alJ 
2009b). The position of LS I +61°303 was found to be 
RA. = 02 h 40 m 22?3, Dec. = 61°13'30"(J2000) with a 95% 
error of 0.069°; in agreement with the accepted position 
dDhawan. M iodusz ewski. & Rupenll2006l) . 

2.1. Spectral Analysis 

The gtlike likelihood fitting tool was used to perform 
the spectral analysis, with "Pass 6 v3" (P6_V3) instrument re- 
sponse functions (IRFs); the P6_V3 IRFs are a post-launch 
update to address gamma-ray detection inefficiencies that are 
correlated with trigger rate. The 10 degree region around 

59 See the FSSC website for details of the Science Tools: 
http://fermi.gsfc.nasa.gov/ssc/datayanalysis/ 




FIG. 1 .— The counts map for 100 MeV-300 GeV in (RA.Dec) of a 10° 
region around the LS I +61° 303 location. The exposure varies by less than 
2.5% across the field at a representative energy of 10 GeV. The source is 
bright and fairly isolated, sitting on a background of Galactic and extragalac- 
tic diffuse emission. A fit to the source yields a significance of more than 70 
a. The dashed line indicates the Galactic equator (b =0); the crosses indicate 
the location of LS I +61° 303 (the brighter source) and a faint nearby point 
source. 




10" icr 

v (Hz) 

FIG. 2. — Fitted spectrum of LS I +61°303 to the phase-averaged Fermi 
data. The solid red lines are the ±1ct limits of the Fermi cutoff power law; 
blue (open circle) data points from MAGIC (high state phases 0.5-0.7); black 
(filled circle) data points from VERITAS (high state phases 0.5-0.8). Data 
points in the Fermi range are likelihood fits to photons in those energy bins. 
Note that the data from the different telescopes are not contemporaneous, 
though they do cover multiple orbital periods. 

the source was modeled for Galactic and extragalactic dif- 
fuse emission, and included one nearby point source at (R.A., 
Dec) of (02 h 23 m 12 s , 62°0'0 "), too faint to be found in the 
3-month Bright Source List dAbdo et al. 2009a|). It is impor- 
tant to include this nearby source in the fitting model because 
at low energies the PSF is sufficently wide that despite be- 
ing ~2.2°away the PSF wings extend across the location of 
LS I +61°303 contributing approximately 13% to the flux 
at this position. Simultaneous modelling of this source ac- 
counts for it's contribution to the flux in this region and any 
uncertainty is folded into the statistical error of the flux of 
LS I +61°303 found by the likelihood fitting tool. The 10° 
region was chosen to capture the broad PSF obtained at 100 
MeV. An alternate fitting method using energy-dependent re- 
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gions of interest was used, yielding compatible results that 
were folded into the systematic errors. 

The Galactic diffu se emission was modeled using GAL- 
PROP, described in IStrong. Moskalenko. & Reimerl d2004l) 
and lStrorigl (12007b . updated to include recent H I and CO sur- 
veys, more accurate decomposition into Galactocentric rings, 
and many other improve ments, including some from com- 
parison with LAT data (lAbdo et alj|2009bl) . The GALPROP 
run designation for our model is 54_59varh7S. The diffuse 
sources contribute ^95% of the observed photons shown in 
Fig.ffl 

Initially a simple power law, E r , was fit to the or- 
bital phase-averaged data yielding a photon index of V ~ 
2.42. However, as indicated in Fig. [2] the energy spec- 
trum appears to turn over at energies ~6 GeV. The pos- 
sibility of an exponential cutoff was investigated, in the 
form E~ r exp[-(E / E cuto fi)]. The chance probability to in- 
correctly reject the power law hypothesis was found to be 
1.1 x 10~ 9 . The best fit exponential cutoff returns a test statis- 
tic (iMattox et al.lll996l) significance value of about 4770, or 
roughly 70ct. The photon index is T =2.21 ± 0.04 (stat) ± 
0.06 (syst); the flux above 100 MeV is (0.82 ± 0.03 (stat) ± 
0.07 (syst)) x 10~ 6 ph cm" 2 s" 1 and the cutoff energy is 6.3 ± 
1.1 (stat) ± 0.4(syst) GeV (see below for a discussion of sys- 
tematics). A total of 135,659 photons were found in the 10° 
region. Evaluating the fit parameters, 6467 ± 80 photons were 
observed fromLS I +61° 303 above 100 MeV. Fig.|2]shows the 
best fit cutoff power law model as well as t he fluxes fit per en 
ergy bin and archival data from MAGIC dAlbert et al 



per en- 
1 120091) 



and VERITAS (lAcciari et alJl2008l) . 

A number of effects are expected to contribute to the sys- 
tematic errors. Primarily, these are uncertainties in the ef- 
fective area and energy response of the LAT as well as back- 
ground contamination. These are currently estimated by using 
outlier IRFs that bracket our nominal ones in effective area. 
These are defined by envelopes above and below the P6_V3 
IRFs by linearly connecting differences of (10%, 5%, 20%) 
at log(E/MeV) of (2, 2.75, 4) respectively. Other potential 
sources of systematic effects investigated are: fitting tech- 
nique; cuts applied (zenith angle, minimum and maximum 
energies); and details of the diffuse modeling. The systematic 
errors estimated using the bracketing IRFs were found to be 
greater than these additional effects, hence the bracketing IRF 
results were quoted for the upper limits on the systematics. 

2.2. Timing Analysis 

LAT light curves were extracted using aperture photome- 
try. The LAT point spread function is strongly energy de- 
pendent and, particularly since LS I +61°303 is located in the 
Galactic plane, there is also significant contribution to the flux 
within an aperture from diffuse emission and point sources 
that depends on the aperture size and the energy range used. 
The aperture and energy band employed were independently 
chosen to maximize the signal-to-noise level. The optimum 
aperture radius was found to be approximately 2.4° in the 
energy range 100 MeV-20 GeV. The time resolution of the 
light curve was 11,478 s, equal to twice the Fermi orbital 
period. Exposures were calculated using gtexposure and 
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FIG. 4. — Power spectrum of the light cur ve. The vertical line indicates 
the known orbital period fromGregory (2002), coinciding with a strong peak 
in the spectrum, while the horizontal lines indicate the marked significance 
levels. 

used to determine the count rate in each time bin. In the expo- 
sure calculation, the spectral shape is assumed to be a power- 
law with a photon index of 2.4. The 1-day binned light curve 
is shown in Fig. [3] Contributions from the nearby source 
and Galactic and extragalactic diffuse backgrounds were esti- 
mated based on the spectral fit and subtracted from the light 
curve 

A search was made for periodic mo dulation by calculat- 
ing th e periodogram of the light curve (lLombll 19761: Icargle 
1982). Since the exposure of the time bins was variable, 
the contribution of each time bin to the power spectrum was 
weighted based on its relative exposure. The periodogram of 
the unbinned, unsmoot hed light curve is shown in Fig. |4] The 
vertical line marks the Gregory (2002) orbital period and a 
highly significant peak is detected at this period. The signif- 
icance levels marked are for a "blind" search with 500 inde- 
pendent frequency steps, however, the effects of the tuning 
of the aperture radius and energy range are not taken into ac- 
count. The period and its error from the LAT observations 
were estimated using a Monte-Carlo approach: light curves 
were simulated using the observed LS I +61°303 light curve 
and randomly shuffling the data points within their statistical 
errors, assuming Gaussian statistics. The corresponding peri- 
odogram was then calculated and the location of the peak at 
^26.5 days recorded. From ^250,000 simulations the distri- 
bution of values gives an estimation of the measured orbital 
period and its associated error of 26.6 ± 0.5 days (lop - 

The binned LAT light curve folded on the Gre gory! (120021) 
perio d with zero phase at MJD 43,366.2749 dGregory et al.1 
1 19791) is shown in Fig. [5] The folded light curve shows a large 
modulation amplitude with maximum flux occurring slightly 
after periastron passage. The overall light curve can be fit 
reasonable well by a simple sine wave, yielding a reduced x\ 
of 1.4 for 1682 d.o.f. However, if w e use the known orbital 
period and ephemeris of the system ( Gregory! 120021) to fit a 
sine wave to each of the individual 9 orbits observed then we 
find that the best fit amplitude varies between 6.8±0.9 and 
2.2±0.9 xl0~ 7 ph cm" 2 s" 1 , which suggests some orbit-to- 
orbit variability. 

2.3. Phase resolved spectral analysis 



FIG. 5. — Folded light curve of LS I +61°303 binned in orbital phase, see 
text for details . The d ashed-lines indicate periastron and apastron as given by 
lAragona et all <200<j) 

The possibility of the spectral shape changing across the 
orbit was explored by running gtlike fits for phase-folded 
bins of 0.1 width. The reduced statistics in each phase bin 
result in a cutoff not being statistically required to fit the data 
and so a simple power law model is used. There is no signifi- 
cant dependence of photon index on phase; a fit to a constant 
value returns a reduced xt °f 1-4 f° r 9 d.o.f, consistent with 
no variation. 

3. DISCUSSION AND CONCLUDING REMARKS 

The Fermi data enable for the first time the detection of a 
modulation in GeV gamma rays at the orbital period of a bi- 
nary system. The derived period is in exc ellent agreemen t 
with the radio and optical-based ephemeris iGregorvl d2002l) . 
The COS-B source 2CG 135+01 is now firmly identified as 
the gamma-ray counterpart to LS I +61°303, resolving a 30- 
year long suspicion that the two were associated. With the 
identification originally based on localization only, the detec- 
tion of orbital-modulated very high-energy emission (>100 
GeV, VHE) from LS I +61° 303 by MAGIC and VERITAS 
dAlbert et alj2006ll2009l:lAcciari et alj2008l) had already pro- 
vided very strong support in favour of this association. 

LS I +61° 303 is detected at a mean flux level above 100 
MeV consistent with that seen by EGRET and AGILE. Aver- 
aged over the orbital modulation, the source persists as one of 
the brightest high-energy g amma-ray sources in the sky over 
a timescale of decades (see lAbdo et al.ll2009al Bright Source 
List, in which this source is the 15 th brightest). The folded 
Fermi light curve peaks around phase 0.3, which is compati- 
ble with periastron passage (when the compact object is clos- 
est to the Be star) acc ording to the latest radial velocity studies 
( Aragona et al. 2009). This contrasts with the behavior at very 
high energies where peak flux occurs at phases 0.6-0.7 and de- 
tections are achieved only at phases ranging from 0.5 to 0.8, 
before or at apastron. In X-rays, LS I +61° 303 also appears 
to pea k at phases 0.6-0.7 (Pa redes et al.lll997t lEsposito et alJ 
120071) whereas the radio peak occurs over a wide range of 
phase s depending upon a 4-year super-orbital cycle (Gregory 
l2002h . 

The average Fermi and EGRET spectra have compatible 
power law indices and fluxes taking into account systemat- 
ics, but the Fermi spectrum also shows a cutoff at approxi- 
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mately 6 GeV. There is no evidence for a phase-dependence 
of the spectral shape and hence, the index or cutoff energy. 
VERITAS reports upper limits during the only VHE obser- 
vations that are contemporary with Fermi, covering only part 
of one orbit from p hase to 0.75 (up to 2008 November 9, 
iHolder et alJ l2009r) . The later phases have short exposure 
times. Moreover, the past VHE history of the source shows 
several non-detect ions at phases 0.6-0.7 (lAcciari et alj|2008t 
lAlbert et"aT] 12009b . perhaps due to variability from one or- 
bital cycle to the other. The Fermi light curve displays signs 
of orbit-to-orbit variability superposed on the mean behav- 
ior, with the primary peak always around phase 0.3. Such 
variability could be attributed to changing conditions in the 
Be star wind, affecting the interaction with the pulsar wind 
or relativistic jet. Indeed, optical spectra show evidence 
for ch anges in wind emission with the orbit (Zamanov et al. 
1999). 

The obvious radiative process to invoke in the HE and 
VHE range is inverse Compton scattering of the abun- 
dant stellar photons into gamma rays by a population of 
electrons accelerated in the vicinity of the compact object 
(e.g. in a relativistic jet or in a pulsar wind). Then, 
all else being equal, the peak flux phase is determined by 
where the seed photon density is highest and by geome- 
try; favorable when the high energ y electrons are seen be- 
hind t he star by the observer, e.g. iDub us. Ceru tti. & Henril 
(120081): iKhangulvan. Aharonian. & Bosch-Ramonl (120081) : 
Sierpowska-Bartosik & Torres (2008). Superior conjunction 
is close in phase to periastron passage in LS I +61° 303 
(0per _ '/'sup = 0.07 to 0.17 depending on the orbital solution). 
Hence, having the Fermi flux peak close to periastron is con- 
sistent with inverse Compton emission from electrons located 
close to the compact object. The cutoff in the average spec- 
trum could arise due to radiative losses (because of differ- 
ent accelerating conditions for electrons, because of the mag- 
netic field amplitude in the relativistic jet or the pulsar wind 
along the orbit and/or because of the greater photon density 
at periastron), or due to a varying maximum energy for ac- 
celerated electrons or to pair p roduction on stellar photons for 
gamma rays above ~50 GeV { [Dubus 2006; Sidoli et al.ll2.006t 
Cerutt LDubus. & Henril2 008: Sier powska-Bartosik & Torresl 
2009). In the latter case, cascade emission might also be 
seen in the Fermi range. All these effects introduce phase- 
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dependent spectral changes. Hadronic interactions related to 
crossings o f the Be star's equatorial wind (d i sk) co uld also 
contribute dChernvakova. Neronov. & Walter! 120061) . This 
would provide an independently varying spectral component 
to explain why the HE and VHE emission peak at different 
phases and vary with orbital cycle. The expectation is that 
hadronic interactions would result in two asymmetric peaks in 
the light curve whose amplitude depends upon the intercepted 
matter density during the crossings and occuring at phases a 
priori unrelated to periastron passage but on the orientation of 
the orbit of the compact object relative to the Be star disk. 

Continued monitoring by Fermi combined with dedicated 
campaigns by pointed instruments is needed to better con- 
strain spectral variability and establish the multiwavelength 
connections: how do orbit-to-orbit variations compare in dif- 
ferent energy ranges? Are there separate HE and VHE spec- 
tral components? 
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